clc,clear
% 载入101个坐标点的经纬度数据（其中1和102都是起点/终点）
sj0 = load('sj.txt');
% 矩阵切片并转化为向量
x = sj0(:,[1:2:8]) ; x = x(:); 
y = sj0(:,[2:2:8]) ; y = y(:); 
%
sj = [x, y]; dl = [70,40]; sj = [dl;sj;dl];
% 角度转化成弧度
sj = sj * pi/180 ;

% 构造距离矩阵
d = zeros(102);
for i = 1:101
    for j = i+1 : 102
        % 根据经纬度和地球半径，将坐标矩阵化为距离矩阵
        d(i,j) = 6370 * acos(cos(sj(i,1)) - sj(j,1) * cos(sj(i,2)) * cos(sj(j,2)) + sin(sj(i,2)) * sin(sj(j,2)));
    end
end
% 刚才求的是三角矩阵，而距离矩阵是对称矩阵
d = d+d';
% 随机数发生器
rand('state',sum(clock));

%  巡航路径 ; 长度
path = []; long = inf; % matlab中inf的大小是21204，是一个足够大的数

% 求较好的初始解（蒙特卡洛方法）：进行一千次循环，每次都更换旅行顺序，直到总旅行距离小于21204
for j = 1:1000
    % randperm：对某个长度的行向量进行整数随机赋值
    % path0是节点的旅行顺序，temp是总旅行距离
    path0 = [1,1+randperm(100),102]; temp = 0;
    % 计算总旅行距离
    for i = 1:101
        temp = temp + d(path0(i), path0(i+1));
    end
    % 检验是否较好（标准有点低了？）
    if temp < long
        path = path0; long = temp;
    end
    % 若效果不好，则更换旅行顺序
end

% e是新解是否接受的判断条件，L是循环次数，at是温度下降比率，T是初始温度
e = 0.1^30; L = 20000; at = 0.999; T = 1;

% 退火过程
for k = 1:L
    c = 2 + floor(100 * rand(1,2)); % 产生新解
    c = sort(c); c1 = c(1); c2 = c(2);
    % 计算代价函数增量
    df = d(path(c1-1), path(c2) + d(path(c1),path(c2 + 1)) - d(path(c1-1), path(c1)) - d(path(c2),path(c2+1)));
    
    % 模拟材料状态改变的概率
    if df < 0
        path = [path(1:c1-1), path(c2: -1: c1), path(c2 + 1:102)]; long = long + df;
    elseif exp(-df/T) >= rand
        path = [path(1:c1-1), path(c2:-1:c1), path(c2 + 1:102)]; long = long + df;
    end
    T = T*at;
    % 接受新解的条件
    if T < e
        break;
    end
end

xx = sj(path,1); yy = sj(path,2);
plot(xx,yy,'- *') % 画出路径图
